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A  USERS’/PROGRAMMERS’  MANUAL  FOR  TWAKE 


1.  INTRODUCTION 


The  computer  program  TWAKE  is  a  modified  version  of  the  CMC-3DPNS  program  described 
in  Refs.  [1-4].  The  parent  code  is  available  in  the  NASA/COSMIC  software  library  (Ref.  [8])  and. 
under  separate  contracts,  variations  of  the  original  code  have  been  acquired  by  the  Naval  Research 
Laboratory  and  the  David  Taylor  Research  Center.  At  NRL  the  program  has  been  used  to  simu¬ 
late  a  variety  of  turbulent,  two  and  three-dimensional,  self-propelled  and  drag  wakes  (Refs.  [5-7]). 
In  addition,  several  unpublished  calculations  of  surface  ship  wakes  and  two-dimensional  turbulent  jets 
have  been  completed.  Comparisons  of  the  calculations  with  experimental  data  have  consistently  shown 
good  quantitative  agreement  with  respect  to  axial  evolutions  of  characteristic  (singie-point)  flow  parame¬ 
ters  such  as  maximum  velocity  deficit  (excess)  and  maximum  turbulence  kinetic  energy.  In  addition,  the 
simulations  of  the  experiments  of  Ref.  [9]  reproduced  the  qualitative  behavior  of  major  flow  structures 
in  the  cross-plane. 

In  order  to  complete  these  numerical  simulations,  it  has  on  occasion  been  necessary  to  alter  and 
merge  portions  of  the  Navy-procured  variants  mentioned  above.  From  these  efbrts  has  evolved  a  version 
of  CMC-3DPNS  which  may  be  described  as  being  “ship  wake  specific”  while  retaining  the  flexibility  of 
the  original  program  through  the  use  of  parameters  defined  during  input. 

NRL  has  been  tasked  to  deliver  the  modified  software  to  the  OCNR  (Code  12)  Ship  Wake  Detection 
Program.  This  report  is  intended  to  serve  as  a  programmers’/users’  manual  in  support  of  that  task.  It 
should  be  considered  as  supplemental  to  Refs.  [2-3]  and  this  report  will  assume  advance  familiarity  with 
that  literature.  The  subsequent  sections  will  describe  the  fundamentals  of  using  TWAKE  to  simulate  a 
class  of  self-propelled  wake  flows  near  a  free  surface.  The  physical  approximations  employed  to  derive 
the  conservation  equations  as  well  as  the  mathematical  concepts  involved  in  their  discretization  are 
described  in  Refs.  [1,4]  and  will  not  be  repeated  herein. 


2.  OVERVIEW  OF  THE  COMPUTATIONAL  MODEL 

The  conservation  equations  to  be  considered  are  the  steady,  three-dimensional,  time-averaged  (in  the 
turbulence  sense)  parabolic  Navier-Stokes  equations.  The  effects  of  turbulence  are  described  using  mod¬ 
elled  transport  equations  for  the  turbulence  kinetic  energy  and  the  isotropic  dissipation  function  along 
with  an  anisotropic  closure  for  the  turbulent  stresses.  Since  the  continuity  equation  is  not  parabolic,  ad¬ 
ditional  formulational  steps  are  necessary  to  render  the  equation  set  amenable  to  a  streamwise-marching 
numerical  solution  technique.  As  is  usual  the  divergence  of  the  transverse  momentum  equations  provides 
an  equation  for  the  static  pressure.  The  pressure  field  that  satisfies  this  Poisson  equation  consists  of 
complementary  and  particular  parts  with  the  complementary  solution  satisfying  the  homogeneous  part  of 
the  original  equation.  An  additional  equation  is  provided  by  the  definition  of  a  harmonic  function  which 
is  forced  to  zero  (within  a  defined  tolerance)  as  the  continuity  equation  is  satisfied.  These  equations  and 
the  ordering  analysis  leading  to  their  derivation  may  be  found  in  Ref.  [4] . 

The  computer  program  has  been  designed  to  consider  the  simultaneous  convection  and  difhsion  of 
up  to  30  arbitrary  scalar  field  variables,  <7,- .  The  general  form  of  the  equation  system  is 


rt  \  dq<  dq'  a  9 
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In  these  equations  n.  is  the  primary  flow  (marching)  direction,  the  subscript  /,  2  <  /  <  3.  denotes  the 
coordinate  directions  defining  the  transverse  solution  plane,  and  n ;  is  the  outward  unit  normal  vector. 
In  Eq.  (1)  Si  is  a  source  or  sink  term  for  q,-  and  x,  is  the  diflision  coefficient.  The  constants  an,i  in 
Eq.  (2)  allow  the  specification  of  the  appropriate  boundary  conditions  on  each  boundary  of  the  solution 
domain.  The  coefficients  c,-  and  d,-  as  well  as  parametric  or  functional  representations  of  the  diflision 
coefficients  and  source  terms  are  supplied  by  the  user  as  part  of  the  input  process  and  by  providing 
additional  coding  as  necessary., 

The  turbulent  PNS  equations  and  modelled  equations  described  in  the  first  paragraph  are  a  subset 
of  the  general  equation  system  and  appropriate  specifications  of  the  coefficients,  functionals,  and  source 
terms  have  been  programmed.  Nevertheless,  there  remains  considerable  latitude  for  the  user,  through 
the  use  of  parameters  specified  on  input,  to  further  tailor  the  equation  set  for  specific  applications, 
especially  with  regard  to  the  treatment  of  source  and  diflision  terms.  Several  of  the  original  parameters 
axe  not  described  in  Refs.  [2-3]  and  some  additional  choices  have  been  added  by  this  author.  The  full 
equation  set  in  a  hybrid,  as  coded,  form  is  included  as  Appendix  A  of  this  report  and  should  serve  to 
clarify  and  complete  the  descriptions  in  the  above  references. 

The  numerical  solution  of  the  governing  equation  set  is  accomplished  by  means  of  a  finite-element 
algorithm.  The  particular  algorithm  employed  by  TWAKE  is  derived  in  Chapter  4  of  Ref.  [4]  using 
the  Galerkin- Weighted  Residuals  formulation.  A  brief  outline  of  the  algorithm  is  useful  to  introduce 
necessary  terminology  and  for  reference  in  subsequent  sections  of  this  report. 

To  a  degree  the  finite-element  algorithm  is  an  integral  transformation,  transforming  an  initial  value 
partial  differential  equation  into  a  larger  order  system  of  ordinary  differential  equations.  The  algorithm  is 
established  by  first  sub-dividing  the  flow  domain  into  a  number  of  sub-domains  or  finite  elements.  Each 
element  is  associated  with  a  number  of  discrete  nodes  located  on  (or  within)  the  element  boundary. 
The  spatial  variation  within  the  element  of  any  flow  variable,  qe,  is  assumed  in  terms  of  interpolation 
polynomials,  N ,  and  the  values  (unknown)  of  the  variable  at  the  node  points,  Q,  for  example, 

qe(xux,)  =  {JVt(*t)}T  {Q(xi)}e  .  (4) 

The  curley  braces  denote  column  matrices  whose  order  is  equal  to  the  number  of  nodes  comprising  the 
element  and  the  elements  of  {Nk}  are  polynomial  functions  complete  to  degree  k.  The  numerically 
determined  finite  element  approximation  qh ,  to  the  true  solution  q ,  is  the  summation  of  the  independent 
solutions  qe, 

M 

q(xuxi)  ~  qh(xi,xi)  =  ^^(x!.*,),  (51 
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where  M  equals  the  number  of  finite  elements  in  the  discretization.  Within  the  Galerkin  formulation, 
for  example,  the  algorithm  statement  for  a  two-dimensional  high  Reynolds  number  boundary  layer  form 
of  Eq.  (A-l)  with  homogeneous  boundary  conditions  is 

5,  J  {Nt}l(u$)dr  =  5e  {Uif  J  {Nk}{Nk}{Nkf  dx2{Ul}' 

ft.  ft. 

+  {Wu}f  J  {Nk)-±-{Nk)~{Nk)Tdx2W,}e  +  {U2}Tc  J  {tVt}{iVl}£-{iVt}rdx2{C'1}4 

ft.  ft. 
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In  Eq.  (6)  S,  is  the  standard  finite-element  assembly  operator,  the  ’‘prime"  denotes  ordinary  differenti¬ 
ation  with  respect  to  js j. ,  and  the  Green-Gauss  theorem  has  been  used  to  transform  the  difhsion  term. 
Each  of  the  integrals  are  analytically  evaluable  upon  specification  of  the  interpolation  polynomials  in  Eq. 
(4).  The  elements  of  the  resultant  hyper-matrices  are  themselves  column  or  square  matrices.  Following 
the  notation  developed  in  Ref.  [4],  Eq.  (6)  is  written 


Se 
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The  meaning  of  the  matrix  notation  is  easily  deduced  with  reference  to  Eq.  (6).  For  instance  [,430b  1  ]  in 
the  convective  term  denotes  that  it  is  the  matrix  formed  by  integrating  over  a  one-dimensional  element 
(4),  the  product  of  3  interpolation  polynomial  matrices  (3),  the  first  two  of  which  are  not  differentiated 
with  the  third  being  differentiated  (001).  Matrices  appropriate  for  the  PNS  equations  have  been  coded 
into  TWAKE  for  linear  interpolation  polynomials  only  ( k  =  1  in  Eq.  (4)),  for  both  one-dimensional 
[4...]  elements  and  two-dimensional  triangular  [5...]  elements. 

Application  of  the  finite-element  algorithm  to  each  of  the  equations  in  Appendix  A  results  in  a 
system  of  ordinary  differential  equations  (and  algebraic)  which  are  compactly  written  as. 


[Cii  m'e  +  [Dii  m.  +  iSi}'  =  {oj 


(3) 


where  {5,}e  contains  all  of  the  non-homogeneous  terms  and  {A}  contains  the  combined  effects  of  con¬ 
vection  and  difhsion.  To  solve  these  equations  the  family  of  one-step  implicit  finite-difference  integration 
algorithms  is  used, 

{*Vi  =  mi+ 1  -  {Qi)j  -  **1  (*  {QiYj+i  +  (1  -  *)  {Q.}>)  =  {0} ,  (9) 


where  j  is  the  n  step  index,  An  is  the  step-size  and  9  is  the  implicitness  factor,  with  9  =  ~  yielding  the 
trapezoidal  rule.  Combining  Eqs.  (9)  and  (8)  to  eliminate  the  derivative  yields  the  non-linear  algebraic 
equation  set, 

{^({Qi}J+i)}.=  {0}.  (io) 

The  Newton  iteration  algorithm  is  applied  to  Eq.  (10) 
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where  p  is  the  iteration  index  and  [7]  is  the  Jacobian 
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The  equations  to  be  solved  at  each  iteration  are  Eq.  (11)  written  in  the  form, 

[j  (Wi}J+i)]  {«*}£!  =  -  {r  {mu t) (13) 

where  {<5Q,},  the  iteration  vector,  is  the  dependent  variable.  The  right  side  of  the  above  equation  is  Eq. 
(10)  evaluated  at  the  pth  iteration. 


W+i  =  [Ci]  (mUi  -  {<?■•}>)  +  An  (*{A}p+l  +(1  -  9){Gi}i)  .  (14) 

where 

{Gi}p  =  [A]  {Qi}p  +  (5,}p  .  i  13 1 

Iterating  the  solution  until  {F}  vanishes  to  within  a  defined  tolerance  implies  the  approximate  vanishing 
of  the  iteration  variable  {SQi}  in  Eq.  (13)  and  hence  convergence.  The  matrix  solution  technique  for 
Eq.  (13)  is  by  L-U  decomposition  and  back  substitution. 
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The  code  has  no  '‘hard-coded”  and  linked  control  sequence  until  full  initialization  has  occurred, 
and  alternate  paths  exist  in  execution  as  well.  The  user  is  thus  required  to  construct  a  problem-specific 

command  sequence  which  will 

1)  generate  the  computational  space, 

2)  specify  the  dependent  variable  set, 

3)  insure  that  fluid  is  initialized  properly, 

4)  specify  appropriate  boundary  conditions  for  each  dependent  variable, 

5)  specify  certain  parameters  which  control  the  flow  path  through  the  execution  module  and. 

fii  choose  among  several  output  options  or  otherwise  provide  output  routines. 

This  has  been  done  in  Refs.  [2-3]  for  a  class  of  semi-bounded  ‘urbulent  flows  with  an  irregular  computa¬ 
tional  domain.  The  next  section  and  major  portion  of  the  remainder  of  this  report  will  explain  an  input 
sequence  appropriate  for  the  simulation  of  surface  ship  wake  flows.  Comparing  the  two  should  efficiently 
educate  the  user  in  the  machinations  of  operating  the  code  and  yield  a  versatility  to  create  additional 
data  sets.  Once  a  command  structure  is  developed  for  a  class  of  flow  geometries,  the  majority  of  the 
previous  discussion  in  this  section  as  well  as  other  characteristics  of  the  program  can  remain  entirely 
transparent  to  the  user.  Since  it  is  likely,  however,  that  improved  modelled  equations  or  constitutive  re¬ 
lations  or  more  efficient  solvers  and  iterative  techniques  will  become  available,  an  outline  of  the  relevant 
sections  of  code  which  might  be  modified  will  be  given  in  Section  4. 

3.  DESCRIPTION  OF  COMMAND  DATA  INPUT 

The  primary  requirement  of  the  user  of  the  program  is  to  construct  an  input  sequence  of  commands 
which  will  describe  the  flow  to  be  numerically  simulated.  Such  a  data  sequence  is  included  as  Appendix 
A  of  Ref.  [2]  for  the  PNS  computation  of  wing-body  juncture  flow.  Another  and  quite  different  input 
sequence  is  included  as  Appendix  B  of  this  report  for  the  turbulent  wake  behind  a  self-propeiled  surface 
ship.  Following  a  few  general  remarks  describing  the  command  data  structure,  the  next  sub-section  of 
this  report  will  explain  the  functions  of  the  fundamental  commands  included  in  Appendix  B  and  indicate, 
where  appropriate,  how  the  data  should  be  altered  to  compute  additional  flow  fields  of  the  same  type. 
This  will  be  followed  by  a  short  sub-section  containing  example  modifications  which  either  allow  effects 
not  included  in  the  Appendix  B  commands  or  significantly  alter  the  geometry.  It  will  be  assumed  at 
this  point  that  the  reader  is  familiar  with  the  texts  of  Appendices  A  and  B. 

The  procedures  for  describing  the  flow-field  geometry,  initial  fluid  state  and  other  characteristics 
pertinent  to  specific  flows  via  the  input  sequence  are  quite  flexible.  At  the  core  of  the  input  module 
is  subroutine  BDINPT.  This  subroutine  sequentially  reads  data  card  images  from  logical  unit  5  (LUoj. 
The  individual  card  images  must  have  the  format  (A8.72A1),  the  first  eight  characters  (left-justified)  of 
which  contain  a  pre-programmed  command  name,  while  the  remaining  characters  contain  any  additional 
parameters  required  by  the  specific  command.  The  subroutines  which  read  these  card  images  translate 
the  literal  data  into  integer  or  real  data  as  indicated.  The  commands  can  be  categorized  into  two  general 
groups,  those  which  operate  on  and/or  store  numerical  data,  and  those  which  cause  control  to  transfer  to 
other  ensembles  of  subroutines  which  may  themselves  require  additional  literal  or  numerical  data  from 
LU5  or  other  logical  units.  The  command  card  images  are  terminated  by  “T”  if  the  data  are  numerical 
and  by  “DONE"  if  the  data  are  literal.  Upon  completion  of  an  operation  specified  by  a  command,  control 
returns  to  BDINPT  which  reads  the  next  card  image.  A  list  of  available  commmands  and  their  functions 
can  be  found  in  Appendix  C  of  Ref.  [2]. 

Immediately  upon  execution  of  TWAKE,  control  resides  in  program  MAIN.  The  principal  functions 
of  this  routine  are  to  specify  the  size,  IZSIZE,  of  the  primary  array  in  the  code.  IZ(IZSIZE).  and  to  call 
subroutine  BDINPT.  Currently  IZSIZE  -  700000  which  is  sufficient  to  compute  a  cross-plane  containing 
51x51  nodes.  Additional  code  has  been  added  to  read  file  names  i  from  LU5)  and  open  tiles  associated 
with  logical  units  required  for  customized  input  and  output  subroutines.  On  some  computer  svstems. 
these  logical  units  may  be  designated  by  way  of  the  job  control  language  and  in  that  event  those  particular 
statements  in  program  MAIN  should  be  deactivated.  The  first  several  card  images  in  the  data  deck  of 
Appendix  B  contain  names  for  these  additional  LU's.  It  is  not  necessary  to  delete  these  card  images 
if  MAIN  has  been  altered  to  attach  these  files  via  job  control  statements.  This  introduces  an  important 
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point  regarding  the  command  data  structure.  If  a  card  image  contains  a  character  string  which  is  not 
recognized  as  a  legitimate  command,  then  control  continuously  passes  to  subsequent  images  until  a  viable 
command  is  detected.  Thus  the  entirety  of  the  annotated  data  set  in  Appendix  B  (or  even  this  entire 
report)  can  be  read  into  the  code  verbatim. 

Discussion  of  Appendix  B 

The  first  sensible  statement  to  BDINPT  in  Appendix  B  is  the  command  FENAME  and  this  command 
must  be  the  first.  The  command  calls  subroutine  FENAME  which  initializes  default  entries  in  the  integer 
and  scalar  arrays,  IARRAY(500)  and  RARRAY(SOO),  respectively.  The  specifications  given  these  entries  in 
the  current  version  of  FENAME  are  appropriate  for  the  surface  ship  wake  calculations  of  Refs.  [6-7]  and. 
in  principle,  the  several  subsequent  commands  in  Appendix  B  of  the  IARRAY  and  RARRAY  type  need  not 
be  present.  They  have  been  included  to  remind  the  user  of  those  entries  which  would  most  likely  be 
altered  to  specify  a  different  scaling  to  a  geometrically  similar  flow.  As  is  apparent  from  Appendix  B  a 
command  of  the  type 
IARRAY  2  7  T 

puts  the  value  of  7  into  the  second  location  of  the  array  IARRAY.  The  comments  in  Appendix  B  suffice  to 
explain  the  meaning  and  use  of  the  several  parameters  which  are  explicitly  referenced.  Descriptions  of 
most  of  the  other  defaulted  entries  in  these  arrays  can  be  found  in  Refs.  [2-3].  They  need  to  be  changed 
only  if  the  user  needs  to  significantly  alter  the  flow  type  and  geometry.  Note  the  specification  of  the 
coefficients  of  the  governing  equations  in  the  form  presented  in  Appendix  A.  Note  also  that  the  last 
several  of  these  array  specification  commands,  along  with  the  ICQND  and  EXIT  commands,  are  preceded 
by  the  letter  “C”  and  consequently  they  have  no  meaning  to  BDINPT.  Removing  the  extraneous  letter  and 
left-justifying,  however,  will  cause  those  particular  elements  of  IARRAY  to  be  assigned  non-zero  values 
and  thus  activate  debug  output  statements  (LU6)  contained  in  the  majority  of  the  subroutines  in  the 
code.  Some  entries  are  assigned  values  greater  than  1  because  those  particular  integers  are  decremented 
by  1  each  time  output  is  requested.  The  ICOND  command  causes  a  print  to  LU6  of  the  current  status 
of  all  non-zero  entries  in  IARRAY  and  RARRAY  and  the  EXIT  will  terminate  execution  when  encountered. 
Obviously  these  parameters  can  be  activated  (I  >  0)  and  deactivated  (I  <  0)  at  strategic  locations  in 
the  input  data  sequence.  There  does  exist  a  NAMELIST  option  in  FENAME  whereby  the  several  parameter 
specifications  in  Appendix  B  can  be  can  be  accomplished  by  direct  reference  to  their  Fortran  variable 
names.  The  option  has  not  been  used  in  Appendix  B  since,  for  illustrative  purposes,  it  is  convenient 
to  isolate  the  parameter  specifications  in  groups  according  to  function.  Most  of  the  more  widely  used 
entries  in  RARRAY  and  IARRAY  are  equivalenced  to  local  Fortran  variable  names  in  FENAME.  The  local 
variable  names  have  been  assembled  into  two  NAMELIST  data  strings,  one  for  each  of  the  arrays.  The 
parameter  specification  in  Appendix  B  can  as  well  be  accomplished  by: 

IARRAY  S00  1  T  NAMELIST  on 

FENAME  T 
JtNAMEOl 

NTYPE*7 ,  N2WAKE=1,  HM=3,  KR0V=19,  LC0L=19,  N0DE=400 

REND 

ANAME02 

UINF=6.7S56,  RH0INF=1.93S,  XMUINF=21 . IE-6 ,  REFL=1 . ,  T0=10., 

TD*20 .  ,  HSINITa.Ol,  HMAX=10. ,  DELP=101. 


The  FEDIMN  command  must  occur  after  the  grid-size  parameters  (NM,  LCOL,  KROW,  NODE  in  Ap¬ 
pendix  B  and  others  defaulted  in  FENAME)  are  prescribed.  This  subroutine  activates  the  linear  finite- 
element  versions  of  the  [A...]  and  [f?...]  matrices  (depending  on  NM)  discussed  in  Section  1  and  loads  default 
entries  into  several  arrays.  The  critical  function  of  FEDIMN.  however,  is  to  partition  the  IZ(IZSIZE)  array. 
The  macro-structure  of  this  array  is  discussed  in  Ref.  [3]  and,  due  to  it’s  critical  importance  to  creating 
new  or  altering  existing  code,  Tables  1  and  2  are  included  to  describe  relevant  further  partitioning.  It 
is  sufficient  at  present  to  point  out  that  the  array  always  appears  as  equivalent  to  it’s  real  counterpart. 
RZ(IZSIZE),  and  every  major  column  matrix  in  the  code,  whether  integer  or  real,  is  a  subset  of  this 
array.  Subroutine  FEDIMN  calculates  addresses  (entry  points  in  IZ/RZ1  for  these  matrices  and  groups  of 
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these  matrices  and  stores  these  addresses  in  the  first  200  locations  of  the  IZ  array.  As  an  example,  the 
following  series  of  Fortran  statements  places  the  value  of  the  r 2  space  coordinate  at  computional  node 
I  into  the  local  variable  Z: 

COMMON  /  ARRAYS  /  IZ (700000) 

DIMENSION  RZC700000),  L(200) 

EQUIVALENCE  (  IZ(1),  RZ(l),  LCD  ) 

EQUIVALENCE  (  L(90),  IX2C0R  ) 

Z  =  RZCIX2C0R  +1-1) 

If  the  debug  parameter,  IARRAY(61),  is  active  upon  encountering  FEDIMN,  these  addresses  will  be  printed. 
The  remainder  of  the  input  sequence  is  devoted  to  initializing  (and  sometimes  further  partitioning)  the 

IZ  array. 

The  next  group  of  command  data  which  deserve  comments  in  addition  to  those  given  in  the  Appendix 
B  annotation  is  the  LINK4  sub-group.  These  commands,  terminating  with  DONE,  construct  the  finite- 
element  domain.  Commands  VX2SCL  and  VX1SCL  determine  node  spacing  in  the  vertical  and  lateral 
directions,  respectively,  by  using  geometric  progression.  The  card  image  immediately  following  VX2SCL. 
for  example,  contains  2  strings  of  numbers,  separated  by  a  comma  for  clarity  (not  necessary),  which 
specify  parameters  for  generating  the  vertical  spacing, 

VX2SCL  T 

0.  9  -l.S  1.  ,  9  -3.  1.  T 

In  this  case  the  sets  of  numbers  specify  2  super-elements  in  the  vertical  direction  such  that  the  first 
super-element  begins  at  xi  =  0.0,  contains  9  finite  elements  (10  nodes  over  the  span),  and  spans  that 
axis  until  x^  =  —1.5,  with  node  spacing  determined  by  a  progression  ratio  ( Pr )  of  1.0.  The  second  super- 
element  also  contains  9  elements  and  continues  the  construction  until  X2  =  —3.0  with  Pr  =  1.0.  Note 
that  the  command  assumes  subsequent  super-elements  begin  where  the  prior  super-element  ended  (-1.5 
in  this  case)  and  the  starting  position  is  omitted.  Since  Pr  —  1.0  (equi-spaced  nodes)  this  construction 
could  have  identically  been  specified  by 
0.  18.  -3.  1.  T 

and  the  first  format  has  been  used  to  remind  the  user  of  the  capability.  If  Pr  <  1.0  the  spacing  becomes 
increasingly  denser  as  the  interval  is  spanned,  for  example,  each  of  the  constructions 
0.  8  -l.S  1.2S  ,  12  -3.  1.  T 

-3.  11  -1.  1.2S  .70.  .8  T 

spams  the  same  domain  (originating  at  different  points)  and  will  cause  a  relatively  denser  spacing  near 
X2  =  0.0.  Similar  comments  apply  to  the  VX1SCL  command.  The  call  to  ELEM  constructs  the  finite- 
element  network  and  node-connectivity  table  using  the  beginning  and  ending  nodes  specified  by  the 
MDECRD  command,  nodes  1  amd  19  in  each  direction,  respectively.  A  schematic  of  the  grid  formed  for  this 
geometry  specification  is  given  in  Fig.  1.  Note  that  since  the  reference  length,  IARRAY(43),  has  been 
assigned  the  value  of  1.0,  the  space  coordinates  correspond  directly  to  ft. 


Fig.  1  -  Schematic  of  finite-element  domain 
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Each  of  the  dependent  variables  has  associated  with  it  an  integer  number  in  accord  with  the  table 
given  in  Appendix  B.  The  XBNQ  sub-group  of  commands  specifies  boundary  conditions  to  these  variables 
according  to  their  respective  integer  designations.  The  default  boundary  condition  is  vanishing  gradient 
normal  to  the  boundary  in  question*  In  Appendix  B,  for  instance,  there  is  no  KBNO  command  for  depen¬ 
dent  variable  5  (turbulence  kinetic  energy)  and  thus  either  dk/dxi  =  0.0  (left  and  right)  or  dk/dxn  =  0.0 
(bottom  and  top)  as  required.  The  domain  created  by  the  ELEM  command  above  corresponds  to  the  free 
surface  wake  experiments  described  in  Ref.  [7].  The  top  boundary  is  assumed  to  be  the  mean  free  sur¬ 
face  in  the  “rigid-lid”  approximation,  the  left  boundary  is  the  symmetry  plane  separating  the  port  and 
starboard  propellers,  and  the  bottom  and  right  boundaries  are  assumed  to  be  tree  of  viscous  influence. 
If  the  initial  (or  boundary)  conditions  are  asymmetric  from  starboard  to  port  then  a  full  planar  solution 
is  necessary.  In  this  case,  for  instance,  the  KBNO  3  command  and  it’s  associated  data  card  image 
would  be  deleted  to  permit  the  lateral  velocity,  U3,  to  float  on  all  boundaries.  The  data  for  o  ‘and  ppi 
would  be  replaced  by, 

KBNO  9  T  PHI  BC:  only  tree  surface  is  impermeable 
LEFT  BOTTOM  2  RIGHT  2 

DONE 

permitting  entrainment  at  all  boundary  planes  except  the  free  surface.  The  “2’s"  in  the  above  data  cause 
the  boundary  node  lists  for  the  bottom  and  right  boundaries  to  begin  at  the  second  node  ‘proceeding 
counter-clockwise)  thus  preventing  the  bounding  of  the  same  node  twice. 

The  next  set  of  commands  in  Appendix  B  beginning  with  DESCRIPT  204  and  continuing  through 
VX3ST  are  concerned  with  specifying  type  of  and  format  for  printed  output  on  LU6.  The  discussions 
in  Refs.  [2-3]  sue  adequate  to  explain  the  general  functions  of  these  statements.  Until  the  command 
IOSAVE  is  encountered  the  data  are  merely  prescribing  headers  for  scalar  output  and  designating  which 
scalars  to  write.  The  DESCRIPT  203  command  provides  header  information  for  vector  output.  There 
is  a  one-to-one  correspondence  between  the  listed  headers  and  the  column  vectors  denoted  by  the  code 
numbers  listed  under  the  IOSAVE  command.  The  IQMULT  command  provides  scalar  multipliers  for  these 
14  vectors  and  in  this  case  the  first  11  are  multiplied  by  the  scalar  in  RARRAY(2)  (which  happens  to 
be  1.0),  the  next  by  RARRAY ( 17 1 ) ;  and  the  last  2  again  by  RARRAYC2) .  These  particular  multipliers 
cause  the  output  variables  to  be  printed  in  their  non-dimensional  form  as  described  in  Appendix  A.  The 
particular  pressure  is  multiplied  by  RARRAY (171)  because  the  code  operates  with  the  variable  divided  by 
that  constant.  As  an  example  to  illustrate  a  couple  of  points,  if  the  user  wants  to  view  the  dimensional 
axial  velocity  and  further  wishes  to  record  the  effective  turbulence  difhsion  coefficient  in  a  form  which 
compares  it  to  the  laminar  viscosity  (see  Eqs.  ( A- 1 ) ,  (A-16)  and  (A-17)),  then  the  last  header  card  under 
DESCRIPT  203  might  be  altered  to; 

U3 ’ U3 ’  PP  P  PHI  NUTRB/NU 

the  last  statement  under  IOSAVE  should  be  altjred  to; 

8248  9248  1247  T 

and  the  IOMULT  data  card  should  be  altered  to; 

27  11*2  171  2*2  47,  15*1  T 

These  alterations,  in  addition  to  designating  the  header,  cause  the  non-dimensional  velocity  variable 
defined  in  Appendix  A  to  be  multiplied  by  RARRAY(27)  (uoo  specified  in  a  prior  command  in  Appendix 
B)  before  printing.  The  multiplication  of  i/<  (variable  number  1247)  by  RARRAY (47)  (reference  Reynolds 
number  calculated  during  the  input  sequence)  produces  the  ratio  of  the  turbulent  and  laminar  difhsion 
in  Eq.  (A-l).  As  to  the  construction  of  the  composite  numbers,  the  material  in  Ref.  [2]  on  pages  21  and 
35  is  contradictory  with  that  on  page  21  correct.  All  arrays  can  be  accessed  using  the  correct  procedures 
of  Ref.  [2]  and  appropriate  addresses  from  the  table  on  pages  55-58  of  Ref.  [3]  and  Tables  1  and  2  of 
this  report.  As  is  explained  in  Appendix  B  the  axial  locations  where  LU6  output  will  occur  are  specified 
by  the  VPVSX  and  VX3ST  commands.  Further  comments  regarding  output  options  including  customized 
output  routines  will  be  reserved  for  Section  4. 


The  LINX  commands  form  the  final  sub-group  which  will  complete  initialization  of  the  rluui.  'Com¬ 
mands  of  this  type  invoke  a  call  to  the  named  subroutine  which  in  turn  accesses  other  modules  to  periorm 
operations  according  to  the  argument.  The  command, 

LINK3  4  T 

for  instance,  as  part  of  the  operations  in  subroutine  LINK3  when  the  argument  is  4.  places  -alls  to 
subroutines  DIMES  and  INDEX.  Subroutine  DIMEN  calculates  several  non-dimensional  parameters  for  later 
use  during  execution  as  well  as  initializes  the  turbulence  model  constants  used  in  the  equations  of 
Appendix  A.  The  call  to  INDEX  finishes  the  partitioning  of  the  IZ  array  begun  in  FEDIMN.  It  is  here  '  hat 
several  of  the  more  commonly  used  addresses  (entry  points  in  the  IZ  array)  are  given  special  storage 
locations  so  that  they  do  not  have  to  be  repeatedly  calculated.  For  example  the  space  between  IZ(48)  and 
IZ(49),  allocated  in  FEDIMN,  is  further  partitioned  into  segments  of  length  NYY*N0DE  I  NYY=IARRAY(90)  i 
The  starting  address  of  each  segment  is  stored  in  the  common  block  JADRES  which  in  most  usages  in  the 
code  has  the  form: 

COMMON  /  JADRES  /  JUl  ,  JU2  ,  JU3  ,  JH  ,  JK  ,  JEPS  ,  JPP  ,  JP  ,  JPHI  ,  JEXC21) 

In  the  execution  module  the  current  values  of  the  dependent  variables  can  be  accessed  by  reference  to 
the  appropriate  element  of  JADRES.  For  instance,  the  value  of  the  non-dimensional  m  velocity  at  node 
I  is  at  RZ( JU1+I-1 ) .  Other  addresses,  which  have  a  one-to-one  correspondence  with  terms  or  groups 
of  terms  in  Eqs.  (9)  through  (15).  •  :e  calculated  and  stored  in  the  common  block  DERIV.  It  should  lie 
noted  here  that  those  users  who  wish  to  modify  code  in  the  execution  module  will  find  a  knowledge  of 
the  code  in  subroutines  FEDIMN  and  INDEX  essential  to  that  task.  Descriptions  of  the  primary  elements 
of  DERIV  can  be  found  in  Table  2. 

The  next  commands,  calls  to  subroutines  GEOHFL  and  NODELM,  complete  the  development  of  element 
matrices  and  vectors  that  are  determinable  solely  from  geometric  data.  It  is  in  GEOMFL.  for  example, 
that  the  derivative  matrix  [A3011]  of  Eq.  (7)  or  it’s  2D-element  counterparts  [3301121  and  [330113]  are 
calculated. 

The  remaining  task  to  be  completed  before  entering  the  execution  module  is  to  initialize  the  fluid 
state  at  the  initial  plane  x\  —  T0.  Inspection  of  the  equation  set  in  Appendix  A  shows  that,  at  the 
minimum,  non-zero  and  positive  starting  values  of  ui,  k,  and  t  must  be  supplied  on  all  nodes  of  the 
discretization.  For  swirling  wake  flows  initial  values  of  uj  and  U3  should  also  be  supplied.  The  user 
will  likely  find  it  convenient  to  develop  his  own  input  routines  to  initialize  the  dependent  var  able  set. 
Examples  of  such  specialized  routines  are  the  subroutines  BLSIUS.  BRDSHW.  CRNINP.  JNCINP.  EDGINP. 
VAKPRO,  VPIDATA,  and  DTHSRDC.  These  routines  establish  a  variety  of  two  and  three-dimensional  bounded 
and  free-boundary  flow  fields  from  either  empirical  formulae  or  by  accessing  pre-processed  data  sets. 
The  principal  requirements  in  routines  of  this  type  are  that  the  ARRAYS  and  JADRES  common  blocks 
be  included  and  that  the  dependent  variables  ultimately  be  specified  .according  to  the  node  ordering 
determined  during  grid  generation,  the  ELEM  sub-group  of  commands  discussed  previously.  While  this 
latter  point  is  obvious,  consider  that  an  alternate  coordinate  specification  to  that  used  in  Appendix  B. 
but  one  that  yields  precisely  the  identical  geometry,  is 


VX2SCL 

T 

-3. 

13. 

0.  1. 

T 

VX1SCL 

T 

3. 

18. 

0.  1. 

T 

This  specification  will  reverse  the  node  sequencing  depicted  in  Fig.  1.  If  the  user  is  uncertain  as  to  how 
the  elements  and  nodes  are  numbered,  activation  of  the  debug  parameters  prior  to  the  ELEM  command 
will  generate  a  listing.  In  this  application  (NTYPE  =  7)  the  LINK2  10  command  accesses  subroutine 
OTHSRDC  which  will  read  free-format  data  from  LU19.  The  local  file  designation  for  LU19  is  given  as 
INPUT19  at  the  beginning  of  Appendix  B  and  this  data  should  be  prepared  in  advance.  In  this  case  the 
first  record  in  INPUT19  is  the  characteristic  dissipation  length  scale.  I4.  Each  succeeding  record  consists 
of  4  entries  corresponding  to  the  starting  values  of  the  dependent  variables  (  ui ,  k.  u;,  113  )  at  each  node  of 
Fig.  1  in  sequence.  The  subroutine  then  initializes  the  dissipation  function,  c,  according  to  the  discussion 
of  Ref.  [7].  The  nodal  data  is  an  interpolation  to  the  computational  plane  of  the  experimental  lata 
of  Ref.  [9],  The  final  two  commands,  LINKS  6  and  LINKS  4.  cause  the  initialization  of  the  laminar 
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and  turbulent  difhsion  coefficients  and  the  Reynolds  stresses  according  to  the  formulation  in  Appendix 
A.  Note  the  presence  in  Appendix  B  of  several  commands  preceding  the  transfer  to  execution  (QKNINT) 
which  are  negated  by  the  letter  '‘C” .  Removing  the  extraneous  letter  will  invoke  a  complete  planar  output 
of  the  geometry  and  fluid  mtialization  for  inspection  prior  to  execution. 

At  this  point  initialization  is  complete  and  the  annotated  data  set  of  Appendix  B  should  be  viewed 
as  the  general  structure  necessary  for  computation  of  surface  wake  flows.  The  following  is  a  list  of  the 
items  in  that  data  set  which  should  be  changed  to  simulate  a  different  flow  of  the  same  general  class: 

1) '  change  the  RARRAY  entries  corresponding  to  the  fluid  properties  (UINF,  RHOINF,  XMUIHF). 

2)  change  the  RARRAY  entries  corresponding  to  the  longitudinal  boundaries  of  the  computational  domain 
(TO,  TO), 

3)  if  there  are  to  be  more  than  19  rows  or  19  columns  in  solution,  change  the  IARRAY  entries  for  KROW 
and  LCOL  to  be  at  least  as  large  as  the  number  of  rows  and  columns,  and  change  NODE  to  be  at  least 

KRCJW*LCOL  +  1. 

4)  change  the  card  images  following  the  VX2SCL,  VX1SCL,  and  NDECRD  commands  to  reflect  the  new 
geometry, 

5)  change  the  KBNO  sub-group  of  commands  to  bound  the  variables  consistent  with  the  new  geometry 
specification, 

6)  change  the  card  image  following  the  VX3ST  command  to  designate  the  longitudinal  stations  for 
complete  planar  output,  and 

7)  change  the  data  on  LU19  to  reflect  the  initial  fluid  state  consistent  with  the  new  geometry  or 
otherwise  supply  a  subroutine  similar  to  DTNSRDC. 

Before  leaving  the  discussion  of  command  data  set  construction  it  is  useful  to  illustrate  by  way 
of  examples  more  general  modifications  necessary  to  calculate  flow  fields  that  deviate  from  the  type 
heretofore  considered.  The  following  sub-section  will  briefly  summarize  modifications  appropriate  for  a 
two-dimensional  flow  and  for  a  flow  with  non-zero  dpjdx i.  Studying  the  contrasts  between  the  data 
sets  (including  that  in  Refs.  [2-3])  will  result  in  an  increased  versatility  in  using  the  code  to  calculate 
dissimilar  flow  fields. 

Additional  Data  Sets 

The  hypothetical  situation  to  be  considered  is  a  two-dimensional  version  of  the  Appendix  B  data. 
In  tnis  case  it  is  supposed  that  the  flow  is  invariant  in  the  lateral  (13)  direction  and  that  the  initial  fluid 
state  is  described  by  the  data  on  the  left  column  of  nodes  in  Fig.  1.  If  the  extent  of  the  computational 
domain  is  the  same,  the  only  necessary  change  to  the  IARRAY  and  RARRAY  specifications  (or  NAMELIST) 
is  to  set  NM  =  2.  It  is  more  efficient  to  alter  LCOL  and  NODE  to  better  fit  the  problem  since  FEDIMN  will 
then  allocate  smaller  blocks  when  partitioning  RZ. 

The  IPINT  vector  should  be  altered  to  delete  the  integration  of  113,  e.g., 

IP  1ST  T 

1562793000  ,  1  -2  -2  -5  6*0  ,  10*11  1  T 

For  the  same  span  in  xj,  the  only  necessary  changes  to  the  ELEM  sub-group  of  commands  are  to  delete 
the  VX1SCL  command  (and  it’s  associated  data)  and  alter  the  NDECRD  command  and  data  to  the  single 
statement, 

NDECRD  -1  T  -1  indicates  2D  grid 

The  IBORD  command  and  data  should  be  replaced  by, 

LINK1  2  T  call  FINDBE:  examine  geometry  for  boundary  elements 

The  boundary  conditions  (KBNO)  should  be  altered  to  delete  the  condition  on  113  and  to  remove  the 
references  to  the  RIGHT  boundary  for  pp  and  <9. 

The  remaining  required  change  concerns  reading  the  initial  fluid  state.  Obviously  subroutine 
DTNSRDC  could  be  altered  to  store  the  first  column  of  the  current  INPUT19  .data  only)  into 

the  RZ  array,  or  a  modified  data  file  could  be  pre-processed  to  be  accessed  by  the  existing  subroutine. 
An  alternative  which  illustrates  further  capabilities  of  the  code  is  to  delete  altogether  the  link  to  TBLINP 
and  DTNSRDC  and  add: 


DEPVAR  1  2  T  load  U1 

1 .01225  1.01193  1.00424  1.01459  1.03879  1.02734  1.01261 
1.01261  1.00843  1.00211  1.00053  9*1.0  T 
RARRAY  2  .01  T  sac  RARRAY(2)=.01 
DEPVAR  5  2  T  load  TKE 

.24124  .24795  .30290  .28394  .25587  .14605  .0430 
.01045  11*1.0  T 

RARRAY  2  .001  T  sat  RARRAY(2) = . 001 
DEPVAR  6  2  T  load  EPS 

.20312  .21166  .28578  .25937  .22188  .09568  .01529 
.00183  11*1.0  T 

RARRAY  21.  T  sat  RARRAY(2)  back  to  1.0 

The  DEPVAR  HI  H2  command  loads  the  RZ  array  beginning  at  entry  points  corresponding  to  the  depen¬ 
dent  variable  specified  by  HI  (see  the  table  in  Appendix  B).  In  this  example  each  indicated  variable  will 
take  on  the  19  values  specified  by  their  respective  data  statements,  with  each  datum  multiplied  by  the 
value  in  RARRAY(H2).  Note  that  the  free  format  allows  a  short-hand  notation  for  repeated  data  and  the 
successive  redesignation  of  RARRAY(2)  disposes  of  the  necessity  to  input  the  leading  zeros. 

The  data  set  in  Appendix  B  of  Ref.  [3]  serves  as  an  example  of  procedures  required  to  include 
non-zero  dpjdxi  into  a  simulation.  The  explanation  given  in  that  reference,  however,  is  extremely 
abstruse  concerning  the  complementary  (inviscid)  pressure  determination  and  only  detailed  deciphering 
of  pressure-related  subroutines  along  with  subroutine  DERVBL  (conservation  equation  assembly)  reveals 
how  to  and  in  what  form  to  communicate  the  gradient  information  into  the  execution  module.  The 
fundamental  requirement  is  to  load  the  pressure  data  into  the  RZ  array  at  locations  IZ(10S+I-1)  where 
I  ranges  over  the  number  of  nodes  in  solution.  Subroutine  DERVBL  expects  the  pressure  gradient  data  in 
the  form: 

dP  (  PooU*  \  d  (  p,  \ 


dxi  \Poa  +  \pooUlcJ  3*1  \Pc *>*%,)  ' 

and  the  user  should  note  that  the  coefficient  is  accessible  as  the  inverse  of  RARRAY(l7i).  Examples  of 
routines  which  accomplish  this  formulation  for  specific  flow  geometries  are  subroutines  PRSGRD.  GETPPR 
and  JHCPPR.  The  first  of  these  calculates  the  pressure  term  for  an  internal  (rectangular  duct)  flow  bv 
global  mass  conservation  and  the  latter  two  compute  the  term  by  interpolation  of  prescribed  axial 
variations  of  pe(*i).  With  the  above  as  background  the  following  word  description  will  help  clarify  the 
operations  performed  by  the  several  commands  on  pp.  91-92  of  Ref.  [3].  The  LINK2  23  command 
reads  and  stores  into  temporary  storage  the  data  which  immediately  follows  that  command.  These 
data  represent  14  groups  of  [xp,  Cp(x3)Xi_x  ]  for  the  geometry  of  Ref.  [1]  under  the  assumption  of 
flow  symmetry  about  the  corner-bisector.  During  the  link  to  subroutine  NQDPCP  these  data  are  used  to 
specify  pressure  boundary  conditions  at  each  node  under  the  KBH0  8  command  (p.  90).  A  solution 
to  Eq.  (A-5)  is  then  obtained  for  each  of  the  14  stations,  the  axial  coordinates  of  which  are  now  stored 
sequentially  beginning  at  RZ(IZ(139)).  The  solutions  are  re-assembled  into  sub-groups  of  14  values 
(solutions)  per  node  and  stored  according  to  the  normal  node  sequence  beginning  at  RZ(IZ(140)). 
During  execution  subroutine  HCDPPR  is  called  at  the  initiation  of  each  marching  step  and  the  xP  interval 
bracketing  the  current  Xi  determines  entry  points  to  the  pressure  table  allowing  evaluation  of  Eq.  (161 
by  linear  interpolation.  These  results  are  loaded  sequentially,  beginning  at  RZ(IZ(10S)),  for  transfer  to 
DERVBL. 

As  a  final  example  which  simply  illustrates  details  of  the  above  description,  consider  that  the  zero 
pressure  gradient  calculation  of  Appendix  B  has  been  established  and  one  wishes  to  account  for  the 
resultant  axial  gradients  of  pp  in  a  subsequent  pass  through  the  flow  field.  In  this  subsequent  iteration, 
the  prior  distribution  of  pp  will  be  identified  with  the  current  pr.  It  is  assumed  that  on  the  first  pass 
the  pp  distributions  were  saved  (unformatted)  on  Lull  in  the  form  specified  in  subroutine  PLT0UT 
For  purposes  of  illustration  it  is  further  assumed  that  these  distributions  were  stored  only  at  the  axial 
positions  specified  under  VX3ST  in  Appendix  B.  however,  it  is  in  the  near  wake  that  the  procedure  has 
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more  relevance.  For  these  conditions  the  following  additional  parameter  specifications  should  be  made 

before  the  call  to  FEDIMN: 

I ARRAY  40 S  1  T  (IPRSCAL)  turn  on  pressure  table  look-up  procedures 

I  ARRAY  394  4  T  (HPoi'AB)  4  tables  oi  pressure  coefficients 

IARRAY  371  361  T  (HPVSXT)  361  values  in  each  table 

IARRAY  161  1444  T  (HPVSX)  storage  allocation  4*361  entries 

The  VPVSX  command  in  Appendix  B  is  removed  and  the  following  series  of  statements  should  follow  the 

VX3ST  specification: 

RARRAY  450  2.  T  temporarily  set  RARRAY(4S0)=2. 

READ  11  161  74  t 

SETVAL  1444  74  74  450  T 

LIHK2  23  T  call  CPSTUP  to  create  tables 

LIHK1  11  T  call  HODPPR  to  initialize  p  grad. 

RARRAY  450  0.  500  T  reset  RARRAY(450)=0 . 

The  READ  command  reads  from  LU11  a  vector  of  length  IARRAY (161)  and  stores  it  temporarily  beginning 
at  RZ(IZ(74) ) .  Since  subroutine  CPSTUP  expects  data  in  pressure  coefficient  form,  the  SETVAL  statement 
multiplies  the  vector  by  RARRAY (450)  and  places  the  result  in  the  same  temporary  storage  location.  The 
call  to  CPSTUP  creates  the  tables  and  finally  HODPPR  initializes  p,  and  dpjdx\. 

The  following  section  will  give  a  brief  overview  of  the  macro-structure  of  the  execution  module  as 
it  applies  to  wake  flows. 

4.  OVERVIEW  OF  THE  EXECUTION  MODULE 


The  QKHINT  command  near  the  end  of  the  Appendix  B  data  transfers  control  into  the  main  drive 
loop.  A  flow  chart  for  the  routines  performing  the  principal  operations  during  execution  is  shown  in 
Fig.  2.  The  flow  paths  branching  to  the  sides  of  QKHIHT  pertain  to  output  options  and  will  be  discussed 
at  the  end  of  this  section.  The  following  several  paragraphs  will  summarize  the  major  functions  of 
the  remaining  subroutines.  Tables  1  and  2  contain  additional  information  regarding  the  RZ  array  entry 
points  for  those  users  having  a  need  for  detailed  code  translation. 

The  call  to  IMPLCT  marks  the  initiation  of  a  step  in  the  axial  direction,  corresponding  to  the  j 
index  of  Eq.  (9).  At  this  point  any  necessary  parametric  evaluations  should  be  performed,  such  as 
those  necessary  for  variable  geometry  or  non-zero  pressure  gradient.  These  types  of  subroutines  are 
user  supplied.  Examples  of  routines  which  fill  the  column  matrix  for  dpc/ dx\  are  subroutines  JNCPPR. 
GETPPR,  and  PRSGRD.  The  first  two  of  these  compute  gradients  by  interpolations  of  either  prescribed 
data  or  solutions  of  Eq.  (A-5).  Subroutine  PRSGRD  calculates  the  pressure  gradient  in  internal  flows  from 
global  mass  conservation  considerations. 

The  sequence  from  DERVBL  through  STRF  in  Fig.  2  is  the  iteration  loop  signified  by  index  p  in  Eq. 
(11).  Subroutine  DERVBL  performs  the  finite-element  assembly  indicated  in  Eqs.  (6)  and  (7).  The  user 
will  recognize  in  DERVBL  the  application  of  the  various  RARRAY  coefficients  discussed  in  Appendix  A  and 
the  ultimate  assembly  of  Eqs.  (10)  and  (14)  for  each  parabolic  equation.  The  master  loop  in  subroutine 
IMPSLV  is  over  the  dependent  variable  index,  HP,  equivalent  to  the  index  i  of  the  equations  in  Section  2. 
As  HP  increments  from  1  to  the  total  number  of  marching  equations  (5  in  this  case),  Eqs.  ( 13)  are  solved 
in  the  order  specified  by  the  IPIHT  command.  Solution  is  by  L-U  decomposition  and  back  substitution 
(BAIDSL).  The  output  from  BAHDSL  is  the  iteration  variable,  6Qi,  forming  the  left  side  of  Eq.  (11).  If  for 
any  variable  at  any  solution  node,  n 
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convergence  is  not  considered  achieved  (any  HPCOHV(HP)  greater  than  zero)  and  <3(’+l  is  determined  bv 
Eq.  (11).  Currently  Hu  (RARRAY(14) )  is  set  to  lO--*  in  FENAME.  Following  IMPSLV.  subroutines  SETDIF 
and  RHLDST  update  the  diFhsion  coefficients  and  Reynolds  stresses.  If  the  convergence  test  was  failed 
on  any  variable,  two  successive  passes  through  subroutines  PPRES  and  STRF  accomplish  the  solutions  of 


the  Poisson  equations  for  pp  and  <z>  and  control  returns  to  DERVBL.If  convergence  was  achieved  for  all 
variables  (all  NPCONV(NP)  equal  zero),  control  returns  to  QKNINT  after  'lie  stress  update.  Flow  charts  of 

DERVBL  and  IMPSLV  are  given  in  Ref.  (3J. 

The  references  to  subroutine  FILTER  in  Fig.  2  deserve  special  mention.  Chapters  5  and  6  of  Ref.  [4] 
discuss  the  possibility  of  waves  occuring  in  the  cross-plane  velocity  component  solutions  for  the  linear 
basis  algorithm  (Jfc  =  1  in  Eq.  (4)).  The  mechanism  for  initiating  these  waves  is  the  action  of  the 
continuity  penalty  function  which  for  the  linear  cardinal  basis  enforces  essentially  a  central  difference 
approximation.  If  a  calculation  is  exhibiting  this  behavior,  simple  grid  refinement  can  exacerbate  the 
problem.  The  original  code  contains  options  to  either  filter  the  secondary  flow  velocity  variables  directly, 
filter  the  iteration  variable  for  the  those  variables,  or  apply  no  filter,  with  the  former  of  these  options 
implemented  by  default.  This  user  has  computed  drag-wake  solutions  (Ref.  [6])  using  both  ••on”  possi¬ 
bilities  for  the  original  filter,  and  observed  only  insignificant  perturbations  in  the  entrainment  induced 
u 2,U3  solutions.  Operation  with  the  filter  off  did  indeed  reveal  the  presence  of  small  spurious  waves  in 
the  <j>  solution  and  a  resulting  infection  of  the  cross-plane  velocities.  It  is  emphasized  that  this  observed 
behavior  was  small  and  sensibly  not  unstable  at  least  for  the  initial  conditions  and  longitudinal  extent 
of  the  computations  in  that  study.  Nevertheless  these  observations  along  with  the  comments  in  Ref.  [4] 
give  reason  to  monitor  and  condition  the  solution  if  necessary.  If  there  is  a  non-zero  swirl  component  in 
the  initial  condition,  however,  application  of  the  filter  to  the  secondary  flow  velocities  is  inappropriate 
since  the  result  is  an  immediate  artificial  smearing  of  those  components.  Grid  refinement  will  alleviate 
the  smearing  effect  at  the  expense  of  a  corresponding  increase  in  computer  time.  Consequently  r.he 
current  code  has  a  filter  option  at  subroutine  STRF  which  will  operate  only  on  the  penalty  function.  The 
default  condition  is  an  active  <i>  filter  set  in  FEHAME  by  IARRAY(403)=1.  The  u^,  U3  filters,  controlled  by 
the  I  ARRAY  entries  (397),  (401),  (402)  and  (404)  are  currently  ofl  The  use  of  these  filter  options  can 
easily  be  deduced  by  examining  the  code  in  subroutines  IMPSLV  and  QKHINT. 

Subroutine  QKHIHT  transfers  control  to  the  output  driver  whenever  the  parameter  IARRAY(86) 
(ICPBT)  is  non-zero.  This  parameter  is  set  on  each  pass  through  IMPLCT  and  will  be  non-zero 

1)  at  values  of  the  step  index  SELF  AS  which  are  integer  multiples  of  KKTPAS  (IARRAY(167)), 

2)  at  the  transverse  planes  closest  to  the  positions  defined  by 

*1  _  nRl3 
R3S  1 00 

where  n  is  am  integer  and  the  RARRAY  entries  are  described  in  Appendix  B,  and 

3)  at  planes  closest  to  the  positions  listed  via  the  VX3ST  command  in  Appendix  B. 

If  KPHT  is  non-zero  full  planar  output  (LU6)  of  all  arrays  specified  with  the  IOSAVE  command  will  result. 
The  format  for  this  output  is  flexible  and  will  automatically  be  determined  to  resemble  the  problem 
geometry  constructed  during  the  ELEM  sequence.  At  axial  stations  where  REOUTP  is  accessed,  subroutine 
PLTOUT  writes  the  entire  dependent  variable  set  (and  Reynolds  stresses)  to  LU15  in  a  column  matrix 
form  corresponding  to  the  node  sequencing  from  grid  construction.  The  nodal  data  is  preceded  by  a 
sequential  listing  of  the  X2,x3  coordinates  of  each  node,  which  for  invariant  geometry  are  written  on 
the  first  pass  only.  If  the  additional  parameters  IL1713  (IARRAY(408))  and/or  ILU8  (I ARRAY (407))  are 
specified  non-zero  during  initialization„then  subroutines  PLTSUR  and/or  STNPRNT  will  be  called  at  each 
axial  step.  Subroutine  PLTSUR  writes  to  LU18  the  full  set  of  computed  variables  at  each  node  in  the 
horizontal  plane  x3  =  0  (top  row  of  Fig.  1),  corresponding  to  the  mean  free  surface.  Currently  a  call 
to  STHPRHT  will  initiate  prints  to  each  of  the  logical  units  8,  9.  16  and  17.  These  are  characteristic  data 
of  the  calculation  and  are  for  the  most  part  comprised  of  extrema  of  the  several  dependent  variables  in 
the  local  transverse  plane.  All  potentially  relevant  IZ  array  entry  points  are  included  in  these  routines 
and  they  can  quite  easily  be  tailored  to  special  needs.  Sample  data  from  all  output  units  are  included 
in  Appendix  C. 
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5.  CLOSING  REMARKS 

The  subroutines  comprising  TWAKE  are  grouped  according  to  their  primary  functions  and  listed 
in  Table  3.  These  modules,  the  input  data  files,  and  all  output  files  resulting  from  the  execution  of  the 
Appendix  B  data  have  been  written  to  magnetic  tape  for  delivery  to  OCNR  (Code  (12)).  The  data  set 
has  been  executed  on  several  NRL  computers  with  identical  results  to  4  significant  digits.  The  following 
is  a  tabulation  of  the  execution  times: 


CRAY 
VAX  11/785 
VAX  11/780 
HP  9000 


2.24 

35.17 

41.83 

215.95 
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TABLE  X  -  Some  important  IZ  array  entry  points 


12 

encry 

commonly 
used  name 

reference 
for  node  I 

remarks 

43 

IC200 

reference  by  element 

standard  matrices:  see  FEDIMN 

48 

IYY 

RZ ( IYY+NYY* ( IPLOC- 1 )* NODE+I- 1 ) 

dependent  variables1 

71 

ITEMP1 

RZ(ITEMP1+(NS+1)*N0DE+I-1) 

Reynolds  stresses2 

74 

ITEMP4 

RZ(ITEMP4+I-1) 

RZCNZF+I-l)  in  IMPSLV4 

89 

IX1C0R 

RZ(IX1C0R+I-1) 

nodal  £3  coordinates 

90 

IX2CQR 

RZ(IX2C0R+I-1) 

nodal  £2  coordinates 

91 

IPRESS 

RZ(IPRESS+I-1) 

nodal  values  of  pr 

92 

IAMU 

RZCIAMU+I-1) 

nodal  laminar  viscosity 

103 

IVEL 

RZCIVEL+I-1) 

temporary  storage,  u2 

104 

IW 

RZ(IW+I-1) 

temporary  storage.  U3 

105 

IPRGRD 

RZ(IPRGRD+I-1) 

nodal  dpjdx\  term 

109 

IGE0M1 

reference  by  element3 

natural  coordinate  derivative  (£3) 

110 

IGE0M2 

reference  by  element3 

natural  coordinate  derivative  1x3) 

136 

IADIF 

RZCIADIF+I-1) 

nodal  values  of  ut 

Notes 

1  IPLOC  is  the  location  in  the  IPINT  vector.  In  Appendix  B,  for  example.  IPL0C=2  for  the  turbulence 
kinetic  energy.  HYY  and  MODE  are  entries  90  and  55  in  IARRAY. 

2  HS={1,2,3,4,5,6}  for  UjU^,  ukuk,  uW3,  U3U3} 

3  See  subroutine  RNLDST  for  examples  of  usage. 

4  For  each  variable  HP  in  turn,  HZF  points  to  the  left-hand-side  of  Eq.  (14)  (including  penalty  term 
for  u2,  U3)  upon  entering  BANDSL.  At  exit  from  BANDSL,  HZF  points  to  left-hand-side  of  Eq.  (11). 

TABLE  2  -  Some  important  addresses  in  COMMON  /  DERIV  / 


DERIV 

entry 

commonly 
used  name 

reference 
for  node  I 

remarks 

25 

HYY 

RZ (IIYY+NYY* (IPLOC- 1 ) *N0DE+I-1) 

IIYY  =  IZ(48) 

36 

HIYB 

RZ(NIYB+NP+I-1) 

dep.  var.1:  {Q,}£+1 

37 

HIZB 

RZCNIZB+NP+I-l) 

{Cily-HL  in  Eq.  (14) 

61 

HIYY 

RZCNIYY+NP+I-l) 

same  as  DERIV  (36) 

62 

HIZZ 

RZ(NIZZ+NP+I-1) 

same  as  DERIV (37) 2 

63 

HIRU 

RZCHIRU+I) 

JADRES( 1 ) — i3 

64 

HIRV 

RZ(HIRV-rl) 

JADRES (2) -1 

65 

HIRW 

RZ(NIRW+I) 

JADRES(3)-1 

71 

HIRE 

RZCNIRE+I) 

JADRES (4) -1 

72 

HIRF 

RZCNIRF+I) 

JADRES (S)-l 

74 

HIYYO 

RZCNIYYO+NP+I-l) 

{<?.};  in  Eq.  (14) 

75 

HIZZO 

RZCNIZZO+HP+I-l) 

{Gi}j  in  Eq.  (14) 

1  HP  =  i  in  Eq.  (1)  for  all  usages  in  table.  HP  for  a  particular  dependent  variable  is  determined 
according  to  the  variable's  location  in  the  IPINT  vector. 

2  RZ(HIZZ+NP+2*N0DE+I-1)  contains  the  first  term  on  right  side  of  Eq.  (14). 

3  JADRES(n)  refers  to  element  of  COMMON  /  JADRES  /  . 
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TABLE  3  -  TWAKE  subroutines  sorted  according  to  primary  function 


Input 

OutDUt 

TJtilitv 

Ean.  Sol 

ADDDEL 

CALORD 

GETDAT 

ASSMAT 

BDATA 

COMOC 

LINK1 

ASMSQ 

BDIHPT 

DRVBUG 

LINK2 

BANCHO 

BLSIUS 

FEPLOT 

LINK3 

HANDSET 

3NDSET 

DPSISO 

LINK4 

BANDSL 

BRDSHW 

ICOND 

LINKS 
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APPENDIX  A  -  PROGRAMMED  FORMS  OF  CONSERVATION  EQUATIONS 

The  programmed  conservation  equations  are  non-dimensional,  with  all  variables  having  been  non- 
dimensionalized  using  reference  values  of  length,  viscosity,  density  and  velocity  as  appropriate,  for  ex¬ 
ample,  _ 

.  _  _  IL  ,/,/•  -  "iua  -  P  t-.J L 

*1  =  r  i  U1  —  '  “lU2  —  2  >  P  ~  o  U-  '  k  ~  (l2  ’ 

Li  (ioo  uooaoo  u  x) 

The  overbars  indicate  conventional  time-averaging  and  the  primes  denote  fluctuating  (turbulence)  com¬ 
ponents.  Omitting  the  asterisks  for  convenience  and  assuming  xi,  x2,  x3  correspond  to  the  streamwise. 
vertical,  and  lateral  directions,  respectively,  the  equations  for  an  isothermal,  constant  property,  incom¬ 
pressible  flow  are  programmed  in  the  following  parameterized  form: 

streamwise  momentum 

dtti  (  dui\  d  |\  i  .  „  .  ,  dui 

Ul^-  +  f?3S5^2— +  u3  — J  -  — [(f?e,  — 
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vertical  momentum 
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lateral  momentum 
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continuity-penalty  function 


complementary  pressure 
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d2Pc  _ 

dxl  dx$ 


Ev*> 


ft  I 

P 

ft 

> 

if, 

%\ 


I 

9. 

y 

v 


I 

IVJ> 

vs 

v5 

•'s' 

s 

$ 

•w 


♦5*1* 

vji 

$ 

W 


vl 

a 


Si 


particular  pressure 

<52PP  ,  32pp 


+  lj  +  i?396  [J:  (db(U2“2)  +  ^(li2“3))  +  J^(J^{u'-u'3)  +  d;<u'3U'3))' 

+^394  [5^(U2£  +  U3fe  +  A398Ul£f)  +  ^(U2fe  +  U3fe  +  /?398Ul£7)] =0,  (A'6) 


where,  p  =  pp  +  pc, 
turbulence  kinetic  energy 


Ul|£-  +  R384  (u7£-  +  U3^)  -  ^  [(^  +  (1  "  ^  £)  I!;, 

"^l[(i?eZl+(1'fl2S3)^)  ^]_?,+e 

_  _  3  (k  f  ,  .  dk  ,  ,  dk\\  d  fk  (  ,  ,  dk  ,  ,  dk  \\ 

+R™Ck  [dT2  (j  + U2“3^j  J +  3^  (j  (U2“35^  +  “3“3^;J  j. 

isotropic  dissipation  function 

<9e  /  de  de  \  „  .  de  /  i/t  de  \  de  ( u,  de 

Ul^7  +  R 384  r2^7  +  U3d^J  " (1  ■  *283)  U d^j  +  dT3  W'dTs 
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=  0,  (A-7) 


-c(lvi  +  c(2j  =  o, 


turbulence  production  . 


■p  _  ,  /  /  dttj  /  /  dt<2  ,  dt*3  ^  /  /  /  /  /  \  /  /  /  /  /  \  *^3 

p  ■  Ul“2ar2  ■  1  3d*3  "  “2“3  ^  “ (  2  2  ■  1  l)  d^  ■  (U3“3 "  “lUl)  557’ 


kinematic  turbulent  stresses: 


k  (  f  du i  V  ,  /dui 


U2U2  —  C3^:  —  l?402C'2l^<-  —  2At03V«  ’ 
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(A- 12) 


( A- 13) 


(A- 1-1) 


(A- 15) 


w*v  w  *wyv  wv 


■»■.' v.'-v.-v,  v."^.- V".  yv WHV  v-<  KY  V.  VY  VjV.  <-. 


The  turbulent  kinematic  viscosity  and  the  reference  Reynolds  number  are  given  by. 


k2 

vx  =  C4 — ,  (A- 16/ 

€ 


and, 


Re  i  = 


Uqq  K 
t^oo 


1  A- 17) 


respectively.  The  function  ^  in  the  cross-plane  momentum  equations  is  included  to  illustrate  the  action 
of  the  penalty  (continuity)  constraint  term  in  the  numerical  algorithm  (see  Ref.  [4j ,  chapter  7).  Aside 
from  the  ''Rnnn"  coefficients,  the  model  contains  the  10  constants  {C i,  C2,  C3,  C4,  c rk,  crt.  Cu .  Cf,  . 
Cic 1  Ce}  which  take  on  the  commonly  accepted  values  {0.94,  0.067,  0.56,  0.068.  1.0.  1.3.  1.44,  1.92.  0.12. 
0.09}.  Most  of  these  model  constants  are  specified  in  subroutine  FEJfAME,  however  a  few  are  computed 
in  subroutines  DIMES  and  DERVBL. 

The  subscripts  on  the  Rnnn  coefficients  in  the  above  equations  refer  to  locations  in  the  scalar  array 
RARRAY(500) .  This  array  is  initialized  to  default  values  in  FENAME  and  modified  to  be  problem  specific 
during  the  input  sequence.  The  coefficients  have  obvious  utility  in  specifying  a  problem  class  from  the 
general  convp576ection-diTh8ipmtion  (Eq.  (1)).  Restricting  attention  to  the  PNS  equations,  several 
of  the  coefficients,  particularly  those  modifying  the  convection  and  source  (pressure  gradient)  terms, 
are  extraneous,  and  have  utility  only  for  debugging  or  sensitivity  studies.  Those  coefficients  modifying 
the  eddy  viscosity  difhsion  and  Reynolds  stress-gradient  terms  have  more  utility  since  they  potentiallyy 
provide  tin  additional  control  over  the  degree  of  non-linearity  and  explicitness  in  the  descretized  equations 
by  specifying  more  or  less  of  the  difiision  in  the  source  terms.  It  should  be  noted,  however,  that  in  this 
version  of  the  code,  the  Reynolds  stresses  are  updated  at  every  iteration,  and  as  such  any  representation 
0  <  Rnnn  <  1  will  give  essentially  identical  results  if  only  the  lead  order  terms  are  retained  in  the 
Reynolds  stress  constitutive  equations.  For  instance,  specifying  #28i  =  1  or  0  in  Eq.  (A-l)  yields  the 
same  result  if  only  the  first  terms  in  Eqs.  (A-13)  and  (A-14)  are  used  and  if  the  stresses  are  computed 
at  the  current  iteration.  In  the  original  version  of  the  code  the  following  default  values  of  coefficients  in 
Eqs.  (A-l)  through  (A-8)  were  set: 

#385  —  #348  =  #282  =  #339  =  #340  —  #396  =  #394  =  #384  —  1 , 

#281  =  #346  =  #399  =  #398  =  #283  =  0, 

#353  =  “  I. 


In  the  current  version  of  the  code  the  above  specifications  are  retained  with  the  exceptions  that  #2 =  0 
and  #346  =  1,  thus  all  diTUsion  terms  are  determined  to  the  same  degree  of  approximation.  Calculations 
were  performed  for  a  variety  of  permutations  of  the  these  coefficients  during  the  studies  reported  in 
Refs.  [6,7]  and  differences  in  the  results  due  to  laminar  difhsion  or  the  alternate  forms  of  Eqs.  (.4-2) 
and  (A-3)  were  insignificant. 

The  stress  specification  in  Eqs.  (A- 10)  through  (A- 15)  is  different  from  the  original  specification  by 
the  presence  of  the  second  terms  in  the  equations  for  u} u’n  and  U1U3.  These  were  added  so  as  to  include 
the  full  model  specified  in  Ref.  [4].  All  of  the  stress  coefficients  have  default  values  of  unity. 
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APPENDIX  B  -  COMMAND  INPUT  DATA  FOR  WAKE  SIMULATIONS 


0UTPUT8 

0UTPUT9 

0UTPUT15 

OUTPUT 16 

0UTPUT17 

OUTPUT 18 

IHPUT19 


FENAME 


IARRAY 


I ARRAY 


TWAKE  DATA  :  DTHSRDC  DESTROYER 
INBOARD  PROPELLER  ROTATION 


INITIALIZE  AT  10  FT. 


Call  SUBROUTINE  FENAME  to  initialize  default  IARRAY  and  RARRAY  entries. 
T  SUBROUTINE  SDINPT  calls  SUBROUTINE  FENAME  on  this  command. 


RARRAY 

RARRAY 

RARRAY 


RARRAY 

RARRAY 

RARRAY 


Define  the  type  of  flow  for  simulation. 

2  7  T  sets  IARRAY ( 2) =NTYPE=7 

NTYPE  controls  subsequent  calls  within  SUBROUTINE  TBLINP  and  SUBROUTINE 
NODPPR.  This  value  will  cause  the  fluid  state  to  be  initialized  according  to 
DTHSRDC  experimental  data  that  has  been  pre-processed  and  saved  in  file 
IHPUT19  above.  This  parameter  also  controls  the  flow  path  for  computing 
the  axial  pressure  gradients  (zero  for  this  case).  See  SUBROUTINE  NODPPR. 

41  1  T  sets  IARRAY(41)=N2WAK£=1 

N2WAKE  alters  the  flow  path  in  the  execution  module  such  that  routines 
appropriate  only  to  flows  near  solid  boundaries  are  not  called. 

Define  the  fluid  parameters. 

27  8.7556  T  sets  the  free-stream  velocity  (UINF)  =  6.7556  ft/sec 
10  1.935  T  fluid  density  (RHOINF)  =  1.935  lbf-sec-sec/(ft4) 

38  21.1E-8  T  fluid  viscosity  (XMUINF)  in  lbf-sec/ft-ft 
Define  the  streamwise  extent  of  computation. 

43  1.0  T  set  reference  length  (REFL)  =  1.0  ft 
24  10.0  T  set  initial  X-plane  (TO)  =10.0 
35  20.  T  set  solution  distance  (TD)  =20. 


The  final  station  will  be  TF=T0+TD=30. 


RARRAY 

RARRAY 

RARRAY 

RARRAY 


RARRAY 


304  0.5  T  set  implicitness  factor  (CIMPTH) 

196  0.5  T  set  implicitness  factor  (THETA) 

7  .01  T  set  initial  axial  stepsize  (HSINIT) 

16  10.  T  set  maximum  allowable  stepsize  (RMAX) 

HMAX  is  expected  as  a  percentage  to  be  applied  to  TD.  In  this  case 
HMAX  will  be  computed  HMAX= . 10*TD=2 . 

13  101.  T  set  the  intervals  for  printing  or  storing  data  (DELP) 

This  is  also  input  as  a  percentage  of  TD  and  in  this  case  will  be  computed 
to  be  greater  than  the  sol'n  distance.  It  will  print  at  the  final  station 
TF=30.  Further  print  controls  will  be  implemented  via  the  VPVSX  and 
VX3ST  command  sequence  used  further  along  in  the  data  set.  Note  that  if 
a  value  less  that  100.,  e.g.,  25.,  had  been  used,  output  would  result  at 
. 25*TD ,  . 50*TD  etc. 
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Define  coaificients  for  conservation  equations:  see  Appendix  A. 

Note  difference  in  RARRAY  command  parameters  to  set  entry  =  0.  See  BDINT. 


RARRAY 

385 

1. 

T  UCMULT 

RARRAY 

348 

1. 

T  VCMULT 

RARRAY 

339 

1. 

T  PCFACT 

RARRAY 

346 

1. 

T  VLDMLT 

RARRAY 

340 

1. 

T  PPFACT 

RARRAY 

396 

1. 

T  OSG 

RARRAY 

394 

1. 

T  T2FIX 

RARRAY 

384 

1. 

T  ECMULT 

RARRAY 

353 

-1 

T  GUMULT 

RARRAY 

282 

0. 

500  T  U2STRS 

RARRAY 

281 

0. 

500  T  U1STRS 

RARRAY 

399 

0. 

500  T  TU1U2P 

RARRAY 

398 

0. 

500  T  T2PFIX 

RARRAY 

283  0. 

Define 

500  T  EXSTRS 

shear  stress  a 

RARRAY 

400 

1. 

T  C112 

RARRAY 

401 

1. 

T  C113 

RARRAY 

402 

1. 

T  C222 

RARRAY 

403 

1. 

T  C223 

RARRAY 

404 

1. 

T  C332 

RARRAY 

405 

1. 

T  C333 

RARRAY 

406 

1. 

T  C122 

RARRAY 

407 

1. 

T  C123 

RARRAY 

408 

1. 

T  C132 

RARRAY 

409 

1. 

T  C133 

RARRAY 

410 

1. 

T  C231 

RARRAY 

411 

1. 

T  C232 

Set  up  finite  element  geometry. 

KM*3  denotes  that  a  3D  parabolic  sol'n  is  to  be  obtained.  2D  triangular 
finite  elements  sill  be  constructed. 

IARRAY  52  19  T  set  no.  of  rows  of  computational  nodes  (KROW)  =  19 

I ARRAY  50  19  T  set  no.  of  columns  of  nodes  (KCOL)  =  19 

IARRAY  55  400  T  set  If  ODE 

NODE  should  be  set  slightly  greater  than  KR0W*KC0L  to  allow  for  sufficient 
storage  allocation  in  dynamic  dimensioning  routine  FEDIMN. 

Removing  the  "C"  from  1st  col.  in  next  3  statements  turns  on  debug  prints. 
CIARRAY  61  1,  86  1  T 

CIARRAY  76  300  ,  7  200,  6  3  T 

CIARRAY  122  1  T 

Removing  the  "C"  in  next  2  statements  will  print  status  of  all  non-zero 
IARRAY  and  RARRAY  entries  to  this  point  and  terminate  execution. 


CICOND  T 
CEXIT  T 


FEDIMN 


Call  SUBROUTINE  FEDIMN  to  allocate  storage  for  variable  length  vectors. 
T 


.VAl.VIAI.' 


Specify  sequencing  for  solution  of  equation  system.  The  following  table 
identifies  the  dependent  variables  by  integer  numbers. 

variable  1  streamwise  velocity  (Ul) 

2  vertical  velocity  (U2) 

3  lateral  velocity  (U3) 

4  stagnation  enthalpy  (H) 

5  turbulence  X.  S.  (TXE) 

S  dissipation  fct.  (EPS) 

7  perturbation  pres.  (PP) 

3  complementary  pres.  (PC) 

9  penalty  fct.  (PHI) 

T  fill  the  IPINT  vector,  read  the  next  card  image 


1  5  8  2  3  7  8  9  0  0 


-3  -3  -S  -5  5*0 


10*11  1  T 


The  first  group  of  10  integers  above  specifies  the  order  in  which  the  eqns . 
are  solved  at  each  iteration,  i.e.,  Ul ,TXE ,EPS , U2 ,U3  etc.  Mote  that  the 
stag,  enthalpy  is  not  computed.  Also  the  complementary  pres,  is  not  computed 
in  this  case  but  is  constant.  The  third  group  of  10  integers,  10*11  1 
(this  string  denotes  the  10  integers  beginning  at  1  with  each  successive 
integer  incremented  by  1;  1,2,3  etc)  re-labels  the  governing  equations 
according  to  the  sequencing  specified  in  the  first  group,  i.e.,  eqn(l)  is 
identified  with  Ul,  eqn(2)  with  TKE  etc.  The  second  group  of  10  integers 
specifies  a  staggered  start  for  the  computation.  Eqn(l)  (Ul)  is  solved 
beginning  at  the  first  axial  step.  The  2nd  and  3rd  eqns.  (TXE  and  EPS) 
are  turned  on  at  step  3.  The  4th  and  5th  eqns.  are  turned  on  at  step  5. 
There  is  no  delay  for  the  remaining  eqns. 

Construct  the  finite  element  network 


LIHK2  14  T  call  SUBROUTIHE  DSCRT2  to  set  up  coordinates. 
VX2SCL  T  compute  node  spacing  for  normal  (X2)  direction 
0. 09-1. 51.  9-3.  1.  T 

VX1SCL  T  compute  node  spacing  for  lateral  (X3)  direction 
0.0,  9  +1.5  1.  93.  1.  T 

SDECRD 

1  19,  1  19  T 

ELEM  T  sat  up  element  connections  for  above  geometry. 
DOME  T  leave  DSCRT2. 


1  19 
ELEM 
DOME 


See  Ref.  2  pp.  8-9  and  Ref.  3  pp.  22-23  for  explanation  of  the  above 

command  sequence. 

This  particular  string  will  result  in  the  construction  of  a  3X3  square 
domain  with  the  origin  of  coordinates  at  the  top  left  corner.  The  X3 
(lateral)  axis  is  spanned  by  19  equi-spaced  nodes  (18  elements)  from 
0.0  to  3.0.  The  X2  axis  is  also  spanned  by  19  equi-spaced  nodes  from 
0.0  to  -3.0.  There  results  361  nodes  and  648  elements. 

Establish  boundary  conditions  for  dependent  variables.  These  command  data 
have  a  prescribed  format.  See  Ref.  2  pp.  21-22,  26  and  Ref.  3  pp .  18-19  for 
a  general  description  of  the  boundary  condition  procedures.  For  any  variable 
the  default  be  is  zero  gradient  normal  to  boundary  in  question.  It  is  only 
necessary  to  list  the  boundaries  (for  each  variable  in  turn)  where  this 
should  not  be  the  case.  For  the  flow  situation  considered,  the  left  boundary 
is  a  plane  of  symmetry,  the  top  boundary  is  a  "rigid  lid"  free  surface 
(equivalent  to  symmetry) ,  and  the  bottom  and  right  boundaries  are  m 
the  free  stream. 


The  next  command  is  a  list  of  boundary  nodes.  ihe  positive  direction  is 
counter-clocic-wise ,  e.g.,  -TOP  creates  a  list  of  the  node  numbers  for  the 
top  row  beginning  at  the  left  most  node  (no.  1)  and  continuing  to  the  last 
node  in  that  row  (no.  19). 


KBNO 

-TOP 

KBNO 

-LEFT 

KBNO 

BOTTOM 

KBNO 

BOTTOM 


IBORD 

LEFT  2  BOTTOM  2  RIGHT  2  TOP  2 

DONE 

LINK1  2  T  call  NODELM  to  establish  element-node  connectivity  table. 

List  variables  and  boundaries  where  value  should  be  ma.ntamed  at  initial 
level . 

KBNO  2  T  Fix  var.  2  (U2)  at  U2=0  on  top  bndry. 

-TOP  DONE 

KBNO  3  T  Fix  var.  3  (U3)  at  U3=0  on  left  bndry. 

-LEFT  DONE 

KBNO  9  T  Fix  var.  9  (PHI)  at  PHI=0  at  freestream  bndry s . 

BOTTOM  RIGHT  2  DONE 

KBNO  7  T  Fix  var.  7  (PP)  at  PP=0  at  fresstream  bndrys. 

BOTTOM  RIGHT  2  DONE 

All  variables  on  all  boundaries  not  specifically  referenced  above  will  be 
subject  to  homogeneous  Neumann  bc's.  These  be ' s  allow  entrainment  from 
the  freestream.  Note  that  ?HI=0  allows  non-vanishing  normal  gradient  and 
mass  eflux  or  influx.  Ul,  TKE,  and  EPS  could  as  well  be  prescribed  as 
fixed  at  the  freestream  boundaries  (in  this  instance)  with  no  significant 
difference  in  the  computed  results. 

Set  up  formats  for  LU6  printed  output. 

The  next  two  commands  store  data  for  the  first  and  last  pages  of  the  LU6 
printed  output.  Anything  placed  after  the  command  and  before  "DONE" 
will  appear  on  the  appropriate  page. 

COMTITLE  T  Read  title  for  last  page. 

3D  WAKE  FLOW  -  DTNSRDC :  DESTROYER  IB  ROT. 

DONE 

DESCRIPT  204  T  Descriptive  title  for  first  page. 

3D  WAKE  FLOW  -  DTNSRDC:  DESTROYER 

DONE 

The  next  set  of  commands  designates  which  scalar  data  is  to  be  printed 
at  the  output  stations  and  specifies  titles  and  scale  factors  to  be  applied 
to  the  data.  See  Ref.  2  pp.  12-15. 

Fill  the  header  title  vector  until  "DONE". 

DESCRIPT  332  T  IOPAR  vector 


REFERENCE 

ENGLISH-FT 

ENGLISH-IN 

M-K-S 

C-G-S 

LENGTH . 

.FT . 

.IN . 

.  M .  . 

CM 

VELOCITY . 

.FT/SEC. .  . 

.N.  A . 

.M/S . 

. CM/SEC. 

DENSITY . 

. LBM/FT3 . . . . 

.N.  A . 

.  KG/M3 . 

.G/CC. . . 

TEMPERATURE. . . . 

. RANKINE .  .  .  . 

.N.  A . 

.  KELVIN . 

.  N.  A 

ENTHALPY . 

. BTU/LBM. . . . 

.  N .  A . 

.KJ/KG . 

N.A 

FR02. SPEC. HEAT. 

. 3TU/LBM-R.  . 

.  N  .  A . 

.KJ/KG-K. . . . 

N.  A 

VISCOSITY . 

. LBM/FT-S . . . 

.N.  A . 

. NT-S/M2 . . . . 

.POISE. . . 

LOCAL  PRESSURE. 

.  PSF . 

.PSI . 

. . U EDGE . 

.CPU  TIME. 

. DPDX1 . . . 

.  .ENERGY. . 

.  INT.  VAR.  . 

NWGEOM  H’S 

.  .H21.  .  . 

.  G23 . 

.  Xl/LREF .  . 

. DX1/LREF . 

.EPSILON. . 

.  DX 1M/LREF 

•V 


>  A*  A*A‘.'*"A  A‘A 


REFL  REY.  NO. 

DONE 

Define  location  in  RARRAY  of  scalars  in  accord  with  above  titles. 

IONUMB  T 

999,  5*200,  999,  200  4*43,  200  27  200  27  27,  200 

10  200  10  10,  200  58  200  58  200,  200  97  200  97 
200,  200  30  200  30  200, 

200  38  200  2*38,  999,  36  2*36  63  300,  100  135  120  200  186  188, 

11  12  14  85  47  T 

Daiina  location  in  RARRAY  of  scale  factors  applied  to  above  scalar  data. 
Positive  values  are  multipliers  and  negative  values  are  divisors. 

MPARA  T 

5*2,  22  162  184  163,  222  164  163,  222  170  174, 

222  165  2,  2  -175  2  2  2 ,  2  2  2  176  2 ,  2  2  2  177  178 , 

-351  2  169  2  2,  3*2  108  2,  5*2,  5*2  T 

Designate  titles  for  profiles  of  output  dependent  variables  until 
"DONE".  See  Ref.  2,  page  15. 


DESCRIPT  203  T 
U1  /  UREF 

U2  /  UREF 

U3  /  UREF 

TKE/TKEREF 

EPS/DISSREF 

Ul'Ul' 

U1'U2' 

U1 '  U3 ' 

U2 ’ U2 ' 

U2 '  U3 ' 

U3'U3> 

DONE 

PP 

P 

PHI 

Define  variables  to  be  output  under  above  titles.  Each  entry  is  a  composite 
number  which  is  decoded  by  the  program.  See  Ref.  2,  pp.  17-21  and  35. 

I0SAVE  T  output  vectors 

1248  2248  3248  5248  6248  3271 

4271  5271  6271  7271  8271  7248 

8248  9248  T 

Define  scalar  multipliers  to  be  applied  to  each  vector  in  turn.  In  this 
case  all  vectors  are  multiplied  by  RARRAY(2)  (which  =  1.0)  except  the 
12th  entry  (PP)  which  is  multiplied  by  RARRAY(171). 

I0MULT  T  scalar  multipliers  for  I0SAVE  vectors. 

11*2  171  2*2  14*1  T 

Specify  axial  locations  for  printed  and  stored  output.  The  program  in 
SUBROUTINE  DPSISQ  designates  a  print  at  axial  locations  where  the 
external  pressure  is  specified.  Even  though  the  freestream  pressure 
is  constant  in  this  flow,  the  print  mechanism  will  be  activated  by 
creating  a  pressure  table. 

The  next  command  specifies  a  pressure  table  with  4  entries,  pressure 
equals  ambient  level  at  each  location. 

VPVSX 

4*2116.8  T 

The  next  command  specifies  the  axial  locations  where  the  pressure  is  at 
the  above  levels.  A  print  on  LU6  will  always  occur  at  the  initial  and 
final  stations. 

VX3ST 

10.  16.  22.  30.  T 

Call  SUBROUTINE  DIMEN  to  compute  non-dimensional  parameters. 

LINK3  4  T  DIMEN  and  INDEX 
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Call  SUBROUTINE  GEQMFL  to  generate  the  element  matrices  and  vectors 
f or  the  natural  coordinate  system.  See  Ref.  4,  Chapt .  3. 

LINK1  3  T  GEOMFL 

Call  SUBROUTINE  NODELM  again  to  compute  element  thickness  and  area  from 
data  calculated  in  GEOMFL. 

LINK1  2  T  NODELM 

Call  SUBROUTINE  TBLINP  (NTYPE=7)  to  link  with  SUBROUTINE  DTNSRDC  and 
initialize  U1,U2, U3.TKE.EPS  on  each  node  in  finite-element  domain. 

Data  will  be  required  from  LU19  (file  INPUT19)  specified  above  and 
prepared  in  advance  by  the  user. 

LINK2  10  T  Distribute  3D  wake  vel.  vector  and  TKE.EPS 

Call  subroutines  SETDIF  and  DFCFBL  to  initialize  diffusion  coefficients 
LINKS  8  T  SETDIF-DFCFBL 

Call  SUBROUTINE  RNLDST  to  initialize  Reynolds  stresses. 

LINKS  4  T  RNLDST 

Initialization  is  complete.  Transfer  control  to  main  drive  loop. 

Removing  "C"  from  1st  col.  in  next  4  commands  will  cause  a  print  of 
all  non-zero  IARRAY  and  RARRAY  parameters,  a  node  table  print,  and  a 
listing  of  the  initial  fluid  state.  The  program  will  terminate  at 
"EXIT"  before  entering  execution  sequence. 

CICOND 

CLINK2  5  T  NODES 

CLINK2  6  T  OUTPUT 

CEXIT 

Turn  off  any  debug  switches  which  may  have  been  turned  on.  Note:  they 
may  be  left  on  during  execution  at  the  expense  of  a  great  amount  of 
output  on  LU6. 

IARRAY  81  0,  86  0  T 
IARRAY  76  0  ,  70,  8  0  T 
IARRAY  122  0  T 

QKNIHT  T  Begin  integration  procedures. 

EXIT  T  Terminate  run  at  end  of  QKNINT  call. 

CASE  END 


APPENDIX  C  SAMPLE  DATA  FROM  OUTPUT  UNITS 
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